data<-read.table("HGUPLUS2_sample_QA",sep="\t",head=T)
dim(data)
data[1,]
summary(data)
pdf("sample.pdf")
hist(data[,1],breaks=seq(0,39,0.01),xlim=c(0,2),ylim=c(0,500),xlab="RTPCR",ylab ="counts",main="RTPCR")
hist(data[,1])
dev.off()
getwd()
pdf("sample.pdf")
hist(data[,1],breaks=seq(0,9604,0.1))
dev.off()
pdf("sample.pdf")
hist(data[,1],breaks=seq(0,9604,1))
dev.off()
pdf("sample.pdf")
hist(data[,1],breaks=seq(0,9604,10))
hist(data[,1])
dev.off()
pdf("sample.pdf")
hist(data[,1],breaks=seq(min(data[,1]),max(data[,1]),10))
hist(data[,1],breaks=seq(min(data[,1]),max(data[,1])+10,10))
dev.off()
pdf("sample.pdf")
hist(data[,1],breaks=seq(0,9610,25))
hist(data[,1],breaks=seq(0,9650,25))
dev.off()
pdf("sample.pdf")
hist(data[,1],breaks=seq(0,9650,25),xlim=c(0,9610),ylim=c(0,50),xlab="Counts",ylab ="Number of barcode calls per sample",main="histogram of Barcode calls"))
hist(data[,1],breaks=seq(0,9650,25),xlim=c(0,9610),ylim=c(0,50),xlab="Counts",ylab ="Number of barcode calls per sample",main="histogram of Barcode calls")
dev.off()
pdf("sample.pdf")
hist(data[,1],breaks=seq(0,9650,25),xlim=c(0,10000),ylim=c(0,40),ylab="Counts",xlab ="Number of barcode calls per sample",main="histogram of Barcode calls")
dev.off()
pdf("sample.pdf")
hist(data[,1])
dev.off()
pdf("sample.pdf")
hist(data[,1],breaks=seq(0,9650,25),xlim=c(0,10000),ylim=c(0,40),ylab="Counts",xlab ="Number of barcode calls per sample",main="histogram of Barcode calls")
hist(data[,2],breaks=seq(0,39000,25),xlim=c(0,39000),ylim=c(0,40),ylab="Counts",xlab ="Number of barcode calls per sample",main="histogram of Barcode calls")
dev.off()
pdf("sample.pdf")
hist(data[,1],breaks=seq(0,9650,25),xlim=c(0,10000),ylim=c(0,40),ylab="Counts",xlab ="Number of barcode calls per sample",main="histogram of Barcode calls")
hist(data[,2],breaks=seq(0,39000,25),xlim=c(0,39000),ylim=c(0,40),ylab="Counts",xlab ="Number of SpikedIn calls per sample",main="histogram of SpikedIn calls")
dev.off()
pdf("sample.pdf")
hist(data[,1],breaks=seq(0,39000,25),xlim=c(0,39000),ylim=c(0,40),ylab="Counts",xlab ="Number of barcode calls per sample",main="histogram of Barcode calls")
hist(data[,2],breaks=seq(0,39000,25),xlim=c(0,39000),ylim=c(0,40),ylab="Counts",xlab ="Number of SpikedIn calls per sample",main="histogram of SpikedIn calls")
dev.off()
pdf("sample.pdf")
hist(data[,1],breaks=seq(0,39000,25),xlim=c(0,39000),ylim=c(0,40),ylab="Counts",xlab ="Number of barcode calls per sample",main="histogram of Barcode calls")
hist(data[,2],breaks=seq(0,39000,25),xlim=c(0,39000),ylim=c(0,40),ylab="Counts",xlab ="Number of SpikedIn calls per sample",main="histogram of SpikedIn calls")
hist(data[,3],breaks=seq(0,39000,25),xlim=c(0,39000),ylim=c(0,40),ylab="Counts",xlab ="Number of Present/Absent calls per sample",main="histogram of Present/Absent calls")
dev.off()
pdf("sample.pdf")
hist(data[,4],breaks=seq(0,1.8,0.001))
dev.off()
pdf("sample.pdf")
hist(data[,4],breaks=seq(0,1.8,0.01))
dev.off()
pdf("sample.pdf")
hist(data[,1],breaks=seq(0,39000,25),xlim=c(0,39000),ylim=c(0,40),ylab="Counts",xlab ="Number of Probesets passed barcode QC per sample",main="histogram of Barcode calls")
hist(data[,2],breaks=seq(0,39000,25),xlim=c(0,39000),ylim=c(0,40),ylab="Counts",xlab ="Number of Probesets passed Spike-in QC per sample",main="histogram of Spike-In calls")
hist(data[,3],breaks=seq(0,39000,25),xlim=c(0,39000),ylim=c(0,40),ylab="Counts",xlab ="Number of Probesets passed Present/Absent QC per sample",main="histogram of Present/Absent calls")
hist(data[,4],breaks=seq(0,1.8,0.01),xlim=c(0,2),ylim=c(0,100),ylab="Counts",xlab ="MAD per sample",main="histogram of MAD")
dev.off()
dim(data)
pdf("sample.pdf")
hist(data[,1],breaks=seq(0,39000,25),xlim=c(0,39000),ylim=c(0,40),ylab="Number of Samples",xlab ="Number of Probesets passed barcode QC per sample",main="histogram of Barcode calls")
hist(data[,2],breaks=seq(0,39000,25),xlim=c(0,39000),ylim=c(0,40),ylab="Number of Samples",xlab ="Number of Probesets passed Spike-in QC per sample",main="histogram of Spike-In calls")
hist(data[,3],breaks=seq(0,39000,25),xlim=c(0,39000),ylim=c(0,40),ylab="Number of Samples",xlab ="Number of Probesets passed Present/Absent QC per sample",main="histogram of Present/Absent calls")
hist(data[,4],breaks=seq(0,1.8,0.01),xlim=c(0,2),ylim=c(0,100),ylab="Number of Samples",xlab ="MAD per sample",main="histogram of MAD")
dev.off()
getwd()
pdf("sample.pdf")
hist(data[,1],breaks=seq(0,39000,25),xlim=c(0,39000),ylim=c(0,40),ylab="Number of Samples",xlab ="Barcode calls per sample",main="histogram of Barcode calls")
hist(data[,2],breaks=seq(0,39000,25),xlim=c(0,39000),ylim=c(0,40),ylab="Number of Samples",xlab ="Spike-in calls per sample",main="histogram of Spike-In calls")
hist(data[,3],breaks=seq(0,39000,25),xlim=c(0,39000),ylim=c(0,40),ylab="Number of Samples",xlab ="Present/Absent calls per sample",main="histogram of Present/Absent calls")
hist(data[,4],breaks=seq(0,1.8,0.01),xlim=c(0,2),ylim=c(0,100),ylab="Number of Samples",xlab ="MAD per sample",main="histogram of MAD")
dev.off()
pdf("sample1.pdf")
data[1,]
cor<- cor(log2(data2[,1]), log2(data1[,3]),method="pearson")
cor<- cor(log2(data[,1]), log2(data[,3]),method="pearson")
main <- "Scatterplot on log2 scale: correlation ="
pas<- paste(main,cor,sep=" ")
plot(log2(data[,1], data[,3], main=pas,xlab="Barcode calls count", ylab="Present/Absent calls count")
abline(lm(log2(data[,1])~log2(data[,3])), col="red")
dev.off()
getwd()
pdf("sample1.pdf")
cor<- cor(log2(data[,1]), log2(data[,3]),method="pearson")
dev.off()
pdf("sample1.pdf")
cor<- cor(data[,1], data[,3],method="pearson")
cor
dev.off()
pdf("sample1.pdf")
cor<- cor(data[,1], data[,3],method="pearson")
dev.off()
getwd()
setwd("/data1/bsi/BORA_processing/devel/Yan_GEO/barcode_present_absent")
data1 <- read.table("result_present_absent")
data2 <- read.table("result_barcode")
data3 <- read.table("result_spiked")
data4 <- read.table("new_result_mad")
data5<-read.table("cutoffs.txt",sep="\t",head=T)
setwd("/data1/bsi/BORA_processing/devel/data_HTHGU/NormSumSrc/HGUPLUS2_GeneSumm_batch_V1/bin")
pdf("scatter_log2.pdf")
vec<-data5[,1]+data5[,2]+data5[,3]+data5[,4]
point<-vec
k<-which(vec>2)
point<-replace(point, k, "x")
k<-which(vec<=2)
point<-replace(point, k, ".")
main <- "Scatterplot on log2 scale: correlation ="
vec<-data5[,1]+data5[,2]
color<-vec
k<-which(vec==2)
color<-replace(color, k, "red")
k<-which(vec==1)
color<-replace(color, k, "green")
k<-which(vec==0)
color<-replace(color, k, "black")
cor<- cor(log2(data2[,2]), log2(data1[,2]),method="pearson")
cor<-round(cor, digits = 3)
pas<- paste(main,cor,sep=" ")
plot(log2(data2[,2]), log2(data1[,2]), main=pas,xlab="Barcode calls count", ylab="Present/Absent calls count",col=color,pch=point)
abline(lm(log2(data1[,2])~log2(data2[,2])), col="red")
vec<-data5[,2]+data5[,3]
color<-vec
k<-which(vec==2)
color<-replace(color, k, "red")
k<-which(vec==1)
color<-replace(color, k, "green")
k<-which(vec==0)
color<-replace(color, k, "black")
cor<- cor(log2(data2[,2]), log2(data3[,2]),method="pearson")
cor<-round(cor, digits = 3)
pas<- paste(main,cor,sep=" ")
plot(log2(data2[,2]), log2(data3[,2]), main=pas,xlab="Barcode calls count", ylab="Spikedin Control calls count",col=color,pch=point)
abline(lm(log2(data3[,2])~log2(data2[,2])), col="red")
dev.off()
pdf("scatter_log2.pdf")
vec<-data5[,1]+data5[,2]+data5[,3]+data5[,4]
point<-vec
k<-which(vec>2)
point<-replace(point, k, "x")
k<-which(vec<=2)
point<-replace(point, k, ".")
main <- "Scatterplot on log2 scale: correlation ="
vec<-data5[,1]+data5[,2]
color<-vec
k<-which(vec==2)
color<-replace(color, k, "red")
k<-which(vec==1)
color<-replace(color, k, "green")
k<-which(vec==0)
color<-replace(color, k, "black")
cor<- cor(log2(data2[,2]), log2(data1[,2]),method="pearson")
cor<-round(cor, digits = 3)
pas<- paste(main,cor,sep=" ")
plot(log2(data2[,2]), log2(data1[,2]), main=pas,xlab="Barcode calls count", ylab="Present/Absent calls count",col=color,pch=point)
abline(lm(log2(data1[,2])~log2(data2[,2])), col="red")
vec<-data5[,2]+data5[,3]
color<-vec
k<-which(vec==2)
color<-replace(color, k, "red")
k<-which(vec==1)
color<-replace(color, k, "green")
k<-which(vec==0)
color<-replace(color, k, "black")
cor<- cor(log2(data2[,2]), log2(data3[,2]),method="pearson")
cor<-round(cor, digits = 3)
pas<- paste(main,cor,sep=" ")
plot(log2(data2[,2]), log2(data3[,2]), main=pas,xlab="Barcode calls count", ylab="Spikedin Control calls count",col=color,pch=point)
abline(lm(log2(data3[,2])~log2(data2[,2])), col="red")
vec<-data5[,1]+data5[,3]
color<-vec
k<-which(vec==2)
color<-replace(color, k, "red")
k<-which(vec==1)
color<-replace(color, k, "green")
k<-which(vec==0)
color<-replace(color, k, "black")
cor<- cor(log2(data3[,2]), log2(data1[,2]),method="pearson")
cor<-round(cor, digits = 3)
pas<- paste(main,cor,sep=" ")
plot(log2(data1[,2]), log2(data3[,2]), main=pas,xlab="Present/Absent calls count", ylab="Spikedin Control calls count",col=color,pch=point)
abline(lm(log2(data3[,2])~log2(data1[,2])), col="red")
dev.off()
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
print(d.AD <- data.frame(treatment, outcome, counts))
glm.D93 <- glm(counts ~ outcome + treatment, family=poisson())
anova(glm.D93)
summary(glm.D93)
glm(counts ~ outcome + treatment,family=gaussian(link="log")
)
fit<-glm(counts ~ outcome + treatment,family=gaussian(link="log"))
coefficients(summary(fit))
savehistory("kk.R")
